Deterministic characterization of stochastic genetic circuits 
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For cellular biochemical reaction systems where the numbers of molecules is small, significant 
noise is associated with chemical reaction events. This molecular noise can give rise to behavior 
that is very different from the predictions of deterministic rate equation models. Unfortunately, 
there are few analytic methods for examining the qualitative behavior of stochastic systems. Here 
we describe such a method that extends deterministic analysis to include leading-order corrections 
due to the molecular noise. The method allows the steady-state behavior of the stochastic 
model to be easily computed, facilitates the mapping of stability phase diagrams that include 
stochastic effects and reveals how model parameters affect noise susceptibility, in a manner not 
accessible to numerical simulation. By way of illustration we consider two genetic circuits: a 
bistable positive-feedback loop and a negative-feedback oscillator. We find in the positive feedback 
circuit that translational activation leads to a far more stable system than transcriptional control. 
Conversely, in a negative-feedback loop triggered by a positive-feedback switch, the stochasticity of 
transcriptional control is harnessed to generate reproducible oscillations. 
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Parallel advances in the conceptual understanding of 
gene regulation along with technological advances in 
molecular biology have given rise to the possibility of 
system-level quantitative kinetic measurements of living 
organisms [l[ and synthetic genetic circuit designs 0, [|[ . 
Interpretation of time-series data from complex networks 
and reliable forward-design of gene circuits depend upon 
detailed quantitative mathematical models 0, 0] . These 
models generally take one of two largely exclusive forms 
- either deterministic formulations with reactant concen- 
tration varying continuously in time and governed by a 
system of rate equations, or stochastic formulations that 
explicitly include the discrete and probabilistic change in 
reactant molecule numbers as each subsequent reaction 
occurs [5] . Both approaches have benefits and associated 
limitations. 

The great practical advantage of rate equation models 
is the ease with which the qualitative behavior of the sys- 
tem can be extracted. By focusing upon the long-term 
behavior, the model dynamics are simplified and one is 
able to gain insight into the expected response of the sys- 
tem Q. Rate equation models, however, neglect the fact 
that chemical reaction networks are composed of species 
that evolve on discrete space - jumping from some num- 
ber of molecules to another as each reaction occurs Q. 
The resulting deviation from the deterministic formula- 
tion is called the intrinsic noise in the system (since the 
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fluctuations arise from the reaction dynamics themselves 
and not from some external source) [1, 0. In cellular 
systems with small numbers of reactant molecules, the 
relative magnitude of the intrinsic noise can be large, 
and can give rise to qualitatively different behavior than 
what rate equation models would predict. A system that 
has several possible stable states, for example, may be 
induced to spontaneous transitions between them as a 
result of intrinsic noise [13, HH , leading to a stochastic 
switching of states. In an excitable system, noise may 
cause oscillations to occur in a model that is otherwise 
stable [l2l [l3l [lij . With a given set of physical param- 
eters, it is possible to simulate explicitly the individual 
chemical reaction events, including the effect of intrin- 
sic noise Nevertheless, the design of synthetic cir- 
cuits, or therapeutics aimed at altering an existing net- 
work, require knowledge of the phase diagram, which 
involves a systematic mapping of the parameter space. 
There, stochastic simulation becomes prohibitively time- 
consuming even for reasonably simple genetic circuits in- 
volving 2-3 genes (see below) , and analytical methods are 
needed. 

A number of analytic studies have been done recently 
to model intrinsic noise in genetic circuits, much of it 
built upon the linear noise approximation [l5| and fo- 
cused upon the noise property itself, e.g., 'noise prop- 
agation' through genetic networks [lg, HZ], the equilib- 
rium distribution of fluctuations about multiple steady- 
states [1 8| and constructive effects of noise in signal pro- 
cessing [H, There has been comparatively little 
work, however, aimed at providing tools to study the 
effect of intrinsic noise on the stability of systems where 
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stochastic models exhibit qualitatively different behavior 
from their deterministic counterparts [2(j. Under these 
conditions, the linear noise approximation alone cannot 
predict qualitative changes in the observable dynamics of 
the system, as for example in the case of noise-induced 
oscillations (2l[. Here we present an analytic method, 
which we call the effective stability approximation (ES A) , 
that extends the applicability of existing deterministic 
methods to include stochastic effects. The method is 
an extension of the linear noise approximation, including 
correction of stochasticity to the deterministic equations 
to the order 1/N (where N is the number of molecules 
in the system). It conveniently connects deterministic 
and stochastic descriptions, allowing systematic explo- 
ration of parameter space while at the same time includ- 
ing the essential effect of intrinsic fluctuations. For the 
two model systems examined here, we find the ESA to 
capture reliably the essential features of those systems, 
correctly estimating the effect of intrinsic noise on the 
phase diagrams of systems dominated by as little as a 
few dozen molecules. 

ESA can be applied to generic models of genetic cir- 
cuits, and a brief tutorial is presented in the Methods 
section with the hope that the approach can be used by 
other investigators to include stochastic effects in deter- 
ministic models. The full mathematical details are pre- 
sented in the Supplementary Material. We illustrate the 
power of the method below by considering two examples 
- an autoregulator with positive feedback (an autoactiva- 
tor) [13] and an excitable genetic oscillator linking posi- 
tive and negative feedback loops [HI, [23[ . The behavior 
of both circuits is conveniently visualized by means of 
a phase diagram that cannot be practically constructed 
using numerical simulations if stochastic effects arc to be 
included. Furthermore, the analysis reveals that the sys- 
tem behavior is completely governed by a few dimension- 
less combinations of model parameters - combinations 
that would be very difficult to infer from simulation data 
alone. We hope that our presentation of the ESA method 
will make it accessible to modelers, bioengineers and syn- 
thetic circuit designers for the analysis of various molec- 
ular circuits, while our description of the behaviors of 
the two model systems will provide quantitative-minded 
biologists with a concrete sense of the effect of stochastic- 
ity as well as a succinct means of characterization (e.g., 
a phase diagram with reduced variables) . 



I. RESULTS & DISCUSSION 

A. Autoactivator 

Perhaps the simplest circuit motif able to exhibit mul- 
tiple stable states is the autoactivating positive feedback 
loop (Figure la) (24|. The circuit consists of a single 
gene encoding an activator. Several autoactivator cir- 
cuits have been experimentally characterized, including 
the autoactivation of CI protein by the Prm promoter of 




FIG. 1: (A) A positive-feedback loop capable of maintaining 
two stable states (2^|. (B) An excitable oscillator that ex- 
hibits noise-induced oscillations [l2, 23]. The autoactivator 
triggers the production of a repressor R that provides neg- 
ative feedback control. (The dashed arrows denote lumped 
transcription and translation, the bold solid arrows denote 
activation, the blunt arrow denotes repression and the wavy 
arrows denote degradation.) 



phage A studied by Isaacs et al. [22j, and the autoacti- 
vation of NtrC by the glnAp promoter of E. coli studied 
by Atkinson et al. [23|]. The autoactivator circuit is ex- 
pected to exhibit either a HIGH state characterized by an 
elevated level of protein synthesis, or a LOW state char- 
acterized by a low basal level of production. We simplify 
the model by assuming that the activator binding and 
mRNA turnover are fast compared to the lifetime of the 
protein activator. The effect of the activator is quanti- 
fied by the activation function g (A/Ka, /) where A is 
the activator concentration, Ka is the equilibrium disso- 
ciation constant of the activator and its cognate binding 
site, and / is the maximum fold-activation in the circuit. 
As a particular example, we assume a Hill-form for the 
activation function g (A/Ka, /), 

with cooperative activation (n = 2) The resulting 
model is a single kinetic equation governing the activator 
concentration A(t) [22l. |25|] (Figure la), 

A A 

— = 7 .g(A)-5.A, (2) 

where 7 is the fully activated rate of protein synthesis and 
S is the protein degradation rate (which in prokaryotes 
is often estimated from the growth rate due to growth- 
mediated dilution). 

In the deterministic limit, when the number of reactant 
molecules is very large, we expect Eq.[5]to adequately de- 
scribe the system behavior. Once initial transients have 
died out, the system will approach a steady-state, and A 
reaches its steady-state value A s where the rate of syn- 
thesis and degradation balance, i.e. 7 • g(A s ) = 5 ■ A s . 
The stability of the steady-state is determined by the re- 
sponse of the system to a small perturbation A p , found 
by linearizing Eq. [5] about A s , 

^ = [7 • g' (A s ) -5}-A p ^\-A p . (3) 
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The expression in the square brackets A = [7 • g' (A s ) —5] 
is a constant that depends upon the model parameters. 
If A is positive, the small perturbations will grow in time 
(A s is an unstable state), while if A is negative, the small 
perturbation will decay (^4 S is a stable state). In the 
stable case, the long-term state of the system can be 
thought of as a point located at the bottom of a valley 
(or basin of attraction) - the more negative the constant 
A, the steeper the valley. As the model parameters are 
varied, the valley may become more fiat (A ~ 0) or even 
develop into a mountain (A > 0), resulting in a loss of 
stability. The parameter space is divided into regions 
of different qualitative behavior (as in Figure 2a, black 
curve); the threshold between these domains indicates 
where A has changed sign and is called the phase bound- 
ary. Although the model seems to depend upon a large 
family of parameters (7,(5, Ka, etc.), the stability of the 
deterministic model is actually described by two dimen- 
sionless combinations of these parameters: the ratio of 
the protein concentration with fully activated promoter 
(A = j/S) to the dissociation constant, A /Ka, and the 
fold-activation, /. 

The effective stability approximation (ESA) we pro- 
pose is an approximation that allows the average effect 
of intrinsic noise to be expressed as a positive correction 
to A, 



A' = A + A c 



(A corr > 0) , 



(4) 



(see Eg. [TBI below). The correction reflects an effective 
flattening of the local landscape by stochastic fluctua- 
tions, making it easier for the system to escape from the 
basin of attraction. Adopting this perspective allows the 
analysis used to study the deterministic model to be ex- 
tended to the stochastic model with only minor modi- 
fication. With A' corrected to include the effect of the 
intrinsic noise, the new phase boundaries are drawn to 
coincide with points in parameter space where A' = 0. 

A major source of intrinsic noise in gene reg ulatory 
networks is so-called translational bursting @, Hyf, where 
each mRNA transcript is translated into several peptides 
before the message is degraded, leading to a burst of pro- 
tein synthesis. Typical values of the 'burst size' b can 
vary from close to zero for poorly translated genes (27[, 
up to several dozen [2^, [2^] depending upon the rate of 
translation and the lifetime of the transcript. When in- 
trinsic noise is included in the autoactivator model, and 
the procedure described in detail in Section III-A of the 
Supplementary Material is applied, we find the correction 
to A is X C orr oc Ah/ A 2 where, 
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is a third dimensionless quantity we call the discreteness 
parameter. This parameter captures the average change 
in protein number when a synthesis or degradation event 
occurs, scaled relative to the protein number required to 
initiate activation jV^ = Ka x V ce u, where Ka is the 



activator dissociation constant and V ee u is the cell vol- 
ume. Increasing the discreteness parameter Ah increases 
the magnitude of the discrete change in activator num- 
bers, and therefore increases the relative magnitude of 
the perturbation to the system caused by the intrinsic 
noise. One would expect the circuit to switch more read- 
ily from stable state to stable state as the magnitude 
of the intrinsic noise is increased, thereby reducing the 
average stability of the circuit. On the other hand, as 
the number of activator molecules increases (Na — > 00), 
the discreteness parameter vanishes and the behavior of 
the system is fully described by the deterministic model. 
Thus, the discreteness parameter Ah represents a distil- 
lation of the complicated effect of intrinsic noise on the 
model behavior, captured in a compact expression that 
would be difficult to extract from numerical simulation 
data. 

As shown in Figure 2a, for the autoactivator the pa- 
rameter space is divided into regions of bistability (two 
stable states) and monostability (one stable state). The 
bistability is most easily lost near the phase boundary 
separating the bistable and monostable states. The cir- 
cuit parameters of Isaacs and co-workers [22[ lie close 
to the left-hand tip of the black triangle in Figure 2a 
(/ « 10), and as they observed in their experiments, the 
noise overwhelms bistability in such a system (c.f. Fig- 
ure 2 A of [12]), leading to rapid transitions between the 
stable states. A much greater fold-activation is required 
to maintain two distinct stable states (as likewise noted 
by the authors). 

Actually, once noise is allowed in the autoactivator 
model, one no longer has stability in the strictest sense 
because there is always a chance that a perturbation will 
switch the system from one steady-state to the other. 
With noise, it is not a question of stability, but rather 
the average escape time from the steady-state [T3, fTl| . 
The longer the escape time (compared with other time 
scales in the problem), the more 'stable' the system. To 
emphasize the effect of the intrinsic noise on the stability 
phase plot, we consider a system with a small number of 
activator proteins (Ka • V C ell = 25 molecules) . Using the 
parameters 7 = 2 protein min -1 , <5 _1 =30 min (a half- 
life of ~ 20 min), Ka = 25 nM and a burst size of b = 10, 
the discreteness parameter in E. coli (V ce u ~ 1/im 3 ) is 
Ah ps 0.2. From Figure 2a (dark gray curve), a max- 
imum fold-activation of / > 40 is necessary to ensure 
long-lived bistable states (shown as a cross on the plot). 
It is possible to explicitly compute the average escape 
time from the stable states for this simple model (see 
[30L |3]| and Section II of the Supplementary Material) . 
Figure 2b compares the average escape time as a function 
of Aq/Ka and / for the case above, with A b = 0.2 (dark 
gray curve in Figure 2a) . Along the dark gray curve, the 
escape time is r = (3±0.5) h, which is about six times 
longer than the protein lifetime (which sets the basic time 
scale of the system's 'memory'). 

The escape time is an indirect measure of the system's 
stability. We have developed a more direct method that 
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FIG. 2: Stability phase plot for the autoactivator (Figure la), including the effect of intrinsic noise. (A) The black dashed 
curve is the phase boundary of the deterministic model with transcriptional activation (Aq/Ka is the fully activated protein 
concentration scaled by the activator/DNA dissociation constant). Increasing the level of intrinsic noise by increasing the 
discreteness parameter Af, (i.e. increasing the 'burstiness' of translation or decreasing the number of molecules) diminishes the 
parameter regime of reliable bistability (Re[\'] < 0). Here, A(, = 0.1 (black solid), 0.2 (dark gray) and 0.3 (light gray). (B) 
The average escape time from the stable state is an indicator of the permanence of the bistability. Here, the dark gray curve 
from Figure 2a corresponds to an escape time of about r = 6, where time has been scaled relative to the protein lifetime <5 _1 . 
(C) As in Figure 2a, but now with translational activation. The range of bistability is considerably widened as transitions from 
the LOW to the HIGH state are supressed. Here, Ka ■ V C ell = 25 molecules and the fully activated burst size is b — 4 (black), 
b — 9 (dark gray) and b = 14 (light gray). 



measures the effective rate of divergence of an ensemble of 
stochastic trajectories. This method is of general applica- 
bility and allows a direct evaluation of the accuracy of the 
ESA. The details of that calculation are reserved for the 
Supplementary Material (see Section III-A.2). Compar- 
ing A' to the effective rate of divergence in the stochastic 
simulations of the autoactivator, the ESA is found to be 
accurate for systems with A& < 0.25. 

The burst size b can be reduced by decreasing the rate 
of translation and indeed Ozbudak et al. suggest that 
many poorly translated genes in E. coli could be the 
result of evolutionary selection against burst noise [27| . 
Alternatively, the method of control in the circuit can 
be shifted from transcriptional to translational activa- 
tion. Although the simple deterministic model remains 
unchanged for either choice of trancriptional or trans- 
lational control, the resulting stochastic model exhibits 
improved stability for translational activation. 

Figure 2c shows the result of putting the translation 
rate under control of the activator. Decreasing the trans- 
lation rate in the LOW state has the effect of shift- 
ing the upper branch of the phase boundary, indicating 
a decrease in transitions from the LOW to the HIGH 
state. The translational autoactivator can tolerate a 
larger range of transcription rates (i.e. higher 7) and 
a lower maximum fold-activation (/ > 20), even for 
large burst size. As above, with 7 = 2 protein min -1 , 
5- 1 = 30 min, K A = 25 nM, V ceU = 1 ^m 3 and a fully 
activated burst size b = 10, a fold-activation of / > 25 is 
required to sustain the bistability (shown as a cross on 
the plot), almost half that required in the transcriptional 
autoactivator above. 

To generate the phase plot for a given stochastic model 
requires division of the parameter space of the model 



(Eq. [2]) into a fine grid, with stochastic simulation per- 
formed at each point. Even after several such simula- 
tions are generated, it is unlikely that the discreteness 
parameter will suggest itself as a key measure of the 
magnitude of intrinsic noise. The ESA method provides 
not only a rapid overview of the parameter space, but 
provides compact expressions characterizing the effect of 
intrinsic noise on the observable dynamics. In the next 
section, we shall apply ESA to the analysis of a more 
elaborate circuit model. 



B. Genetic oscillator 

Oscillating systems underlie many physiological pro- 
cesses in the cell, from circadian rhythms [32] to the cell 
cycle itself 33]. In addition to the natural systems, sev- 
eral synthetic genetic oscillator designs have been stud- 
ied, including the mutually-repressing ring-oscillator (Re- 
pressilator) of Elowitz and Leibler [34| and the activator- 
repressor design of Atkinson and co-workers [23[ (which 
has a great deal in common with the model discussed 
below). A recurring motif in experimentally character- 
ized networks is a negative feedback loop serving as a 
system reset [lH, [35|. Without some time delay or in- 
tervening mechanism to prevent reversibility, the system 
will rapidly approach an intermediate equilibrium, and it 
is found both theoretically [3t| and experimentally [33j 
that a negative feedback loop alone is not sufficient to 
maintain reliable oscillations. If, however, the feedback 
repressor is controlled by a bistable autoactivator, the 
oscillations become more robust and coherent since the 
bistable switch acts as a ratchet that 'locks' into the 
HIGH state generating a large amount of repressor to 
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feed back and reset the system to the LOW state where 
the system remains until the activator accumulates over 
a critical threshold to initiate another cycle [ji?}- This 
motif is highly represented in natural gene networks [35| , 
and we shall use the ESA to ascertain the contribution of 
intrinsic noise to the performance of such an oscillator. 

We consider the generic model proposed by Vilar 
and co-workers to describe circadian rhythms in eukary- 
otes [l2|, with a transcriptional autoactivator driving 
expression of a repressor that provides negative con- 
trol by sequestering activator proteins through dimeriza- 
tion [13L 1 3 81 ] . The repressor and activator form an inert 
complex until the activator degrades, recycling repres- 
sor back into the system. In their model, the degrada- 
tion rate of the activator, 5a, is the same irrespective of 
whether it is bound in the inert complex or free in solu- 
tion. We simplify their original model somewhat, and as 
in the previous section, we assume fast activator/DNA 
binding and rapid mRNA turnover, leading to a reduced 
set of rate equations governing the concentration of acti- 
vator A, repressor R and the inert dimer C, 



dA 
~dt 
dR 
~dt 
dC 
~dt 



• \ ■ 'i ( j - 1 • .1- '<'(• • * ' !>' 

( -^•/« - <>/,■ ■ /?- *<■■ A - ll + <) 1 • C 



k c ■ A - R- 5 A ■ C. 



(6) 



We further assume no cooperativity in activator binding 
(n = 1 in the activation function g) and the nominal 
parameter set used in [!§]■ For this more complicated 
system, there is a larger number of dimensionless com- 
binations of parameters that characterize the system dy- 
namics. The scaled repressor degradation rate e = Sr/Sa 
is a key control parameter in the model since oscillations 
occur in the deterministic system only for an intermedi- 
ate range of this parameter. For the nominal parameter 
set used in the deterministic model exhibits oscil- 

lations over the range 0.12 < e < 40 (Figure 3a, black 
region). We shall focus on the parameter regime near to 
the phase boundary at e rs 0.12 and examine the role in- 
trinsic noise plays in generating regular oscillations from 
a deterministically stable system. 

Applying the ESA to the oscillator model, the parame- 
ter Ab A = (pA + t) / (2- K A'Vceii) emerges as an important 
measure quantifying the discreteness in activator synthe- 
sis (see Eq. 36 in the Supplementary Material). Here 
again, bA is the burst size in the activator synthesis, Ka 
is the activator/DNA dissociation constant and V ce u is 
the cell volume. (Here, V ce ii = 100/im 3 as is appropriate 
for eukaryotic cells.) 

Using the nominal parameter set of Vilar et al. [12j in 
our reduced model leads to a burtiness in activator syn- 
thesis of bA — 5 (giving = 6 x 10~ 2 ) and a burstiness 
in repressor synthesis of bn = 10. The phase boundary 
predicted by the ESA is shown as a solid line in Figure 
3a, bounding a region of parameter space between the 



deterministic phase boundary where qualitatively differ- 
ent behavior is expected from the stochastic model. We 
examine the system behavior in this region by running a 
stochastic simulation using the parameter choice e = 0.1 
and Ab A = 6 x 10~ 2 (denoted by a cross in Figure 3a). 
With this choice, the deterministic model is stable (Fig- 
ure 3b, black line). Nevertheless, a stochastic simula- 
tion of the same model, including protein bursting and 
stochastic dimerization, clearly shows oscillations (Figure 
3b, dotted line). 

The time between successive peaks in the stochastic 
simulation of Figure 3b is denoted by T. As is clear from 
Figure 3b, T is itself a random variable. Each simula- 
tion run generates a collection of inter-spike times from 
which the mean (T) and the variance (((T) — T) 2 ) 1 / 2 
can be calculated. Following Steuer et al. [131 ]. the qual- 
ity of the noise-induced oscillations is measured using 
the noise-to-signal ratio rj T = (((T) - T) 2 ) 1 / 2 / (T) , and 
the system is said to exhibit regular oscillations where 
r/T is small [TH, [38J. The dependence of t]t on the re- 
pressor degradation rate e is shown in Figure 3c, with 
the discreteness parameter Af, A = 6 x 10~ 2 (as in Fig- 
ure 3b), using at least 200 spikes to calculate t\t- At 
low repressor degradation rate, the noise-to-signal ratio 
is high, indicating large variance in the inter-spike time T 
and corresponding to a stable {i.e., non-oscillatory) sys- 
tem. As the repressor degradation rate is increased, the 
variance in the inter-spike time T decreases with a con- 
sequent decrease in the noise-to-signal ratio i]t, indica- 
tive of a more regularly oscillating system. Physically, 
the intrinsic noise in this parameter range is sufficient to 
drive the system away from the deterministically stable 
steady-state, yet the noise is not so strong that the return 
trajectory through phase space is much affected. 

As in the autoactivator model, it is useful to compare 
the phase boundary predicted by the ESA to some in- 
dependent measure of stability, in this case rjr- In Fig- 
ure 3c, the ESA phase boundary (for A^ A = 6 x 10~ 2 ) 
is denoted by the interface between the white and gray 
regions, corresponding to a value of tjt ~ 0.2. Using 
data such as that shown in Figure 3c, the points in the 
phase plot with ryr — 0.2 can be found for a range of dis- 
creteness parameter A& A (Figure 3a, filled circles) . These 
points correspond very well to the phase boundary calcu- 
lated using the ESA (Figure 3a, solid line). The results 
are as one would expect - near the deterministic phase 
boundary, very little molecular noise is required to sus- 
tain oscillations, and reasonable periodicity persists even 
for small values of the discreteness parameter {A\, A — > 0, 
bu 7^ 0) . As the repressor degradation rate e is decreased 
to a region favoring stability, more noise is required to 
overcome the deterministic stability of the system and 
initiate the autoactivator trigger. It is illustrative to re- 
mark that each data point in Figure 3a, obtained from 
stochastic simulation 5] , took roughly a day to generate 
on a dual processor desktop computer since at low repres- 
sor degradation rate, a large separation of timescales is 
introduced necessitating long stochastic simulation runs 
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FIG. 3: (A) Stability phase plot as a function of the scaled repressor degradation rate e = 5r/5a for the circuit shown in Figure 
lb. The discreteness in the activator synthesis, Ab A , characterizes the average discrete change in activator concentration during 
each reaction, and consequently the magnitude of the intrinsic noise. The intrinsic noise expands the region of instability (gray) 
extending the parameter range over which oscillations are expected to occur. The deterministic phase boundary is located at 
e » 0.12 (dashed line separating the black and gray regions). The solid line is the phase boundary predicted from the roots of 
Eq. [12] and filled circles denote the phase boundary found by stochastic simulation (see text). The model and parameters are 
as in Vilar et al. 12]. (B) The circuit exhibits noise-induced oscillations (dotted line) with inter-spike time T. The parameters 
used in the simulation correspond to a deterministically stable system (black line). Numerical simulation data was generated 
using Gillespie's direct method Q, with parameters as used in [12| and e = 0.1, Af, A = 6 x 10~ 2 (cross in Figure 3b). (See 
Section III-C of the Supplementary Material.) (C) A plot of the noise-to-signal ratio r\T = {((T) — T) 2 ) 1 ^ 2 / (T) as a function 
of e. The oscillations are regular when tjt is small (the region of noise-induced oscillations predicted by the ESA is gray), and 
rjT was calculated using at least 200 spikes for each point. 



to capture the slowly- varying dynamics of the system. 
By contrast, the solid line generated from the roots of 
Eq. took less than an hour to produce on the same 
machine. Thus, even for a two-gene circuit with several 
degrees of freedom, the ESA affords a compact and con- 
venient means to survey the phase space, drawing atten- 
tion to those regions of particular interest that may be 
probed in more detail by more realistic (though also more 
computationally costly) stochastic simulation methods. 

II. METHODS 

The effective stability approximation can be applied 
to generic models of genetic circuits in a straightforward 
way. Here, a brief outline of the method is provided. 
A self-contained tutorial on stochastic modeling and the 
ESA is found in the Supplementary Material. 

A useful abstraction of genetic regulatory networks is 
as a system of ordinary differential equations [3^, [40| . 
(Here, and throughout, we shall assume a spatially ho- 
mogeneous environment.) We denote the concentrations 
of the reactants of interest by the state vector x, where 
the Xi correspond to the concentration of mRNA, tran- 
scription factors, protein products, etc. The kinetic equa- 
tion governing the evolution of the system takes the form 
^ = f(x), where f is a vector of nonlinear functions 
of the state variables. We can estimate the long-time, 
or steady-state, behavior of the model by first comput- 
ing the equilibrium points x s that satisfy the algebraic 
constraint f(x s ) = 0. We then Taylor expand the reac- 
tion rate vector about the equilibrium point by making 



the substitution x = x s + x p (where x p is an infinitesi- 
mal perturbation away from x s ), and retain only linear 
terms in x p . The resulting dynamics of x p are given 
by^Xp = J • x p , where J is the Jacobian or response 
matrix: = dfi/dxj. The eigenvalues of J are the 
matrix analogue of the parameter A introduced in Eq. [31 
and in a similar fashion if the eigenvalues all have neg- 
ative real-part, then x s is a stable steady state. (There 
are, of course, limitations to how far one can trust the 
linearization [6[, but for our purposes it is sufficient as a 
first approximation.) 

To include stochastic effects in the mathematical 
model, chemical reaction rates must be re-written in 
terms of the reaction propensity and stoichiometry [f|. 
For example in the positive autoactivator example above, 
with the individual synthesis and degradation stoi- 
chiometrics written explicitly, the deterministic model 
equations (Eq. read, 

bursty synthesis: A -A* A + b; v\ = ^ • g (A), 
linear degradation: A ^> A — 1; v% = 5 ■ A. 

We encode this information concisely as the propensity 
vector v = [vi, v%[ — [7 • g(A)/b, 5 ■ A] and the stoichiom- 
etry matrix S = [6,-1]. The discrete change in molecule 
numbers following the completion of a chemical reaction 
causes a deviation from the deterministic solution since 
the deterministic model assumes an infinitesimally small 
and continuous change in the state. (Consequently, the 
deterministic model only applies to systems with large 
numbers of molecules.) We denote the deviation of the 
stochastic model from the deterministic model by the 
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fluctuating quantity to ■ a(t), where u> = l/\/V ce u and 
a(t) describes the stochastic deviation in each species x. 
The \/V ce u scaling arises from the observation that the 
relative magnitude of the intrinsic noise scales roughly as 
the inverse square-root of the number of molecules [la ] . 
Elf and Ehrenberg [2lJ have developed an algorithmic 
expression for the statistics of a using the linear noise 
approximation of van Kampen [l5l |. In that formulation, 
the mean and covariance of the fluctuations about the 
deterministic state are written compactly in terms of the 
propensity vector v and the stoichiometry matrix S; here, 
we shall apply their method to characterize the fluctua- 
tions about the stable state. The first step in the cal- 
culation of the moments of the fluctuations Loa.(t) is to 
construct the auxiliary matrices T and D, evaluated at 
the stable state x s , 

Tijit) = ^ir i = ^r D = s - dia sM- sT (8) 

The drift matrix T = J is the response matrix (or Ja- 
cobian) described above and reflects the local stability 
of the deterministic system to small perturbations [411 ]. 
The diffusion matrix D captures the strength of the fluc- 
tuations and is related to the magnitude of the reaction 
step-size [2l|, [42| ■ It is straightforward to show that to 
leading-order in to the mean of the fluctuations is zero 
((a) = 0) and the variance, denoted by the symmetric 
matrix S = (a ■ a T ), is determined by the solution of the 
system of algebraic equations [l5[ , r • H + E • T T + D = 0. 
Since the fluctuations about the stable state are station- 
ary, the time autocorrelation function depends upon the 
time difference only, and is given by the matrix exponen- 
tial, 



(a. (t) ol t (t - t)) = exp [Tt] ■ 3. 



(9) 



The effect of the fluctuations on the deterministic steady- 
state is calculated by including an additional term in the 
deterministic linearization above: x = x s +Xp+wa. Lin- 
earizing J in oj, we have a stochastic differential equation 
governing the decay of the perturbation modes x p , 



— x - rj<°) 



J (1) W1'X P 



(10) 



The fluctuations affect the decay of the infinitesimal dis- 
turbance well as the dynamics of the average (x p ), 
which (provided u> 3^'(t) <C J^ ^) is approximately gov- 
erned by the convolution equation [43|, [44| , 



- <x p (t)) = JW (x p (t)) + u; 2 Jj c (t- t) (x p (r)) dr, 

o 

(11) 

where J c (t-r) = /jW (t) e J (0) (t-r)j(i) ( r ) J) is ma d e 

up of linear combinations of the cross-correlations 
(ai(t)aj(r)) given by the i row and the j th column of 



the right-hand side of Eq.[9j In the noiseless case, the sta- 
bility of the perturbation x p is determined by the eigen- 
values of diag{Ai} = P 1 • J^ ^ • P where the ma- 
trix P is made of the eigenvectors of 3^ . The analogues 
of the eigenvalues for the convolution equation above are 
found from the poles of the Laplace transform, denoted 
A', which solve the resolvent equation [ia] . 



det 



A'l - J(°) 



1 



V C ell 



Jc(A') 



0, 



(12) 



Here to 2 has been replaced by V c Jl and J c (s) = 

oo 

J 3 c (t)e~ st dt is the Laplace transform of J c (t). If the 

o 

deterministic eigenvalues are distinct, we can further ap- 
proximate the effective eigenvalue A^ by, 



A' = A 4 + 



1 



[P- 1 - J c (Ai)-P]i 



(13) 



where 



denotes the i diagonal entry of the matrix. 



Physically, we interpret the leading-order noise correc- 
tion as the power in the fluctuations at eigenfrequency 
Ai projected in the eigendirection of A;. Since the cor- 
rection term is quadratic, it is always positive and thus 
(ie-stabilizes the eigenmode upon which it is projected. 
(Hence, in Eq. 4 we write A corr > 0.) 

It often happens that out of the term 1/V ce u [ P _1 • 
J c (Xi) ■ P ]a there appears a small parameter that quan- 
tifies the effect of the intrinsic noise. (For the two exam- 
ples above, the small parameters are A;, and Ah A , each 
characterizing the discreteness of the protein change.) In 
the limit that this parameter goes to zero, the effect of 
the intrinsic noise becomes negligible, at least in that 
particular eigenmode. 

Finally, the ESA can be easily implemented in a sym- 
bolic computational environment, without attending to 
the mathematical details (see Section IV of the Supple- 
mentary Material) . A version of the ESA coded in Math- 
ematica is freely available from the authors by request. 
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III. SUPPLEMENTARY MATERIAL 



about the equilibrium point: x = x s + x„ 



Much theoretical work has been devoted to quantifying 
the conditions under which microscopic fluctuations have 
macroscopic effects [46]. The most useful results are of- 
ten restricted to systems with a single degree of freedom 
or employ sophisticated tools such as Ito's calculus. In 
what follows, we aim to develop a convenient and simple 
scheme to assess the stability properties of a dynami- 
cal system subject to molecular noise described by the 
chemical Master equation. The method is an extension 
of the familiar linear stability analysis of nonlinear dy- 
namical systems, although here the effective eigenvalues 
about the equilibrium points are adjusted to reflect the 
influence of the noise. 



IV. MATHEMATICAL METHODS 

A very useful qualitative picture of the behavior of a 
system of nonlinear differential equations emerges from 
the linearized dynamics about the fixed-point(s) (also 
called the steady-state(s)) of the system, defined as the 
reactant concentrations at which the synthesis and degra- 
dation rates balance. The stability of the system near the 
fixed-points can be estimated by calculating the eigenval- 
ues {Xi} of the resulting linearization, which are gener- 
ally a set of complex numbers. If the real parts are all 
negative, we say the system is locally stable, meaning 
small perturbations away from the steady-state are au- 
tomatically corrected. 

Since genetic circuits, both natural and engineered, 
rely upon transfer of information through small num- 
bers of molecules, significant fluctuation is simply one of 
the inherent operating conditions [47| . resulting in noise 
that may give rise to behavior that is very different from 
the behavior predicted by deterministic models. Conse- 
quently, for cell-scale modeling we propose to modify the 
deterministic notion of stability by calculating the effec- 
tive eigenvalues A^, which include the averaged influence 
of the intrinsic noise, 



dt 



J(0) 



(15) 



A' — Xi + X c 



(14) 



Here X corr cx V~„ is inversely proportional to the cell 
volume V ce u For notational convenience in the following, 
we introduce a parameter to that is related to the cell vol- 
ume by: ui~ 2 — Vceii- Sometimes cj~ 2 is called the 'sys- 
tem size', expressing as it does the relationship between 
reactant concentration and molecule numbers [9, MM ■ 



A. Stochastic stability equation 



dt 



To calculate the stability of the macroscopic model 



(Here, and henceforth, we adopt the convention of writ- 
ing all matrix variables in bold upper-case, and all vec- 
tors in bold lower-case.) The eigenvalues of the Jacobian 
j( ' = _ provide the decay rate of the exponen- 
tial eigenmodes; if all the eigenvalues have negative real 
part, we say the system is locally asymptotically stable. 
We shall restrict ourselves to this notion of stability, al- 
though it does ignore algebraically gro wing modes which 
may be important in certain instances [48j . 

To accommodate fluctuations on top of the small per- 
turbation x p , we set x = x s + x p + Lua(t). The Jacobian 



<9x 



will then be a (generally) nonlinear function of the fluc- 
tuations about the steady-state at(t). (As a technical 
aside, we note that we are justified in replacing x by 
x s + x p + uja(t) in both the right- and left-hand side of 
the deterministic model 4| = f(x) since the fluctuations 
at(t) have non-zero correlation time (as we show below) 
and zero mean, allowing us first to conclude that the 
time-derivative of a(t) exists and further that the aver- 
age of this derivative must vanish: (t?) = ^jp- = 0). In 
the limit u> — > 0, we can further linearize J with respect 
to u, 



83 



The stability equation is then given by, 



|x p = [j(°)+o;jW(t)]-x p . 



(16) 



This is a linear stochastic differential equation with ran- 
dom coefficient matrix 3^\t) composed of a linear com- 
bination of the steady-state fluctuations ot(t) which have 
non-zero correlation time (c/. Eq. [9]). We therefore need 
not appeal to any specialized calculi (e.g. Ito's calcu- 
lus) for interpretation since the non- vanishing correlation 
time of the fluctuations ensures that diffcrcntiable 
process and the equation falls under the purview of ordi- 
nary calculus [49l | . 

Our present interest is in the mean stability of the equi- 
librium point. Taking the ensemble average of Eq. I16[ 



iw=' (,) 



(*) 



The right-most term is the cross-correlation between the 
process x p and the coefficient matrix J^ 1 ^). Since the 
correlation time of jW(t) is not small compared with the 
other time scales in the problem, it cannot be replaced 
by white noise, and an approximation scheme must be 



f (x) to small perturbations, the system is linearized developed to find a closed evolution equation for (x p 
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B. Bourret's mode-coupling approximation 

By assumption, the number of molecules is large so the 
parameter u) is small, although not so small that intrinsic 
fluctuations can be ignored. To leading-order in lo, the 
trajectory x p (t) is a random function of time since it is 
described by a differential equation with random coeffi- 
cients. Derivation of the entire probability distribution of 
x p (t) is usually impossible, and we must resort to meth- 
ods of approximation. We shall adopt the closure scheme 
of Bourret [43], 0, to arrive at a deterministic equa- 
tion for the evolution of the averaged process (x p (t)} in 
terms of only the first and second moments of the fluc- 
tuations. In that approximation, provided J*- -* 3> luJ^, 
the dynamics of (x p ) are governed by the convolution 
equation, 



- (x„ (i)) - J (x„ (t)) (17) 
t 

+w 2 J J c (t-r) (x p (r))dr, 
o 

where J c (t - r) = [t) e J<0) (*- T )jM (r)\ is the time 

autocorrelation matrix of the fluctuations and e J ° T is the 
matrix exponential. The equation can be solved formally 
by Laplace transform, 

<x p (s)) = \sl - J(°) - oj 2 3 c (x„ (0)) , 



where now J c (s) = J J c (t) e~ st dt. A necessary and suf- 
o 

ficient condition for asymptotic stability of the averaged 
perturbation modes (x p (t)) is that the roots A' of the 
resolvent, 



det 



A'l - J - oj 2 3 c (A') 



0. 



(18) 



all have negative real parts (Re(X') < 0) [45|, \EM- Some 
insight into the behavior of the system can be gained 
by considering a perturbation expansion of the effective 
eigenvalues A' in terms of the small parameter w. We 
further diagonalize diag[A;] = P" 1 • J^ -* • P, and 

provided the eigenvalues are distinct, we can explicitly 
write AJ in terms of the unperturbed eigenvalues A^ to 
O(oj 4 ) as, 



A' = A, + 



U> 



J c (A l 



(19) 



where [ ■ ]a denotes the i th diagonal entry of the matrix. 

Notice the matrix product 3 c (t — r) contains lin- 
ear combinations of the correlation of the fluctuations 
(ai(t)aj(r)), and as such we must derive an expression 
for those moments. 



Calculating the statistics of the steady-state 
fluctuations 



The statistics of the fluctuations a are fully deter- 
mined by the solution of the chemical Master equation 
(defined below) that comes from treating each reaction 
event probabilistically. In that probabilistic formulation, 
our state at any time t is represented by the vector of 
molecule numbers n € N d ; with m representing the num- 
ber of molecules of a given species. Each reaction causes 
a transition from the initial state n to some new state n' 
reflecting the addition or removal of molecules by that 
reaction. The probability that the transition n — > n' oc- 
curs is the product of the probability of being in state n 
at time t, P(n, f), and the transition probability of mov- 
ing from n — > n', denoted by W n ^ n >. We thus write the 
probability conservation as a balance of flux into and out 
of the state n, which yields a discrete differential equation 
for P(n,t), 



dP (n, t) ^ 



dt 



W n ^„P(n',i)-W n ^P(n,f). (20) 



The evolution equation for P(n, t) is called the Master 
equation [52|. It is rare that the Master equation can be 
solved exactly for P(n,t), and approximation schemes 
are required. One such scheme, the linear noise approx- 
imation llal , is versatile and will be described briefly 
(see also [2 llj and [13]). The approximation begins with 
the assumption that the molecule concentrations can be 
meaningfully separated into a component that evolves 
deterministically, which we shall denote x(i) , and fluctu- 
ations a(i) that account for the deviation of the stochas- 
tic model from the deterministic model. We introduce 
a scaling parameter lj, where uj~ 2 = V ce n is the volume 
of the cell and is an extensive measure of the number 
of molecules. We then make the ansatz that the fluctua- 
tions scale as the square-root of the number of molecules: 
u! 2 ni = Xi + to cti [HI, HH . In that way, a perturbation 
expansion as the number of molecules gets large (u> — > 0, 
with concentration held fixed), returns to zero'th order 
the macroscopic reaction rate equations, 



rfx 



/(*)• 



(21) 



The first-order equation, that comes at 0(u>), charac- 
terizes the probability distribution for the fluctuations 
Tl(a,t) centered on the macroscopic trajectory x(t), and 
has the form of a linear Fokker-Planck equation, 



dU 
~dt 



= -J2 r^(a;ii) + \ ]T /m„ii. 



(22) 



where di denotes d/dai and 

r «(*) = ^ £> = S • diag^] • S T , (23) 
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(see main text). The matrices T and D are independent 
of a, which appears only linearly in the drift term. As 
a consequence, the distribution H(a,t) will be Gaussian 
for all time. In particular, at equilibrium the fluctuations 
are distributed with density, 



'27rrdetH 



exp 



1 



and variance H — (a ■ a T ) determined by, 

r-E + H-r T + D = o. 



(24) 



Furthermore, the steady-state time correlation function 



(a (t) a T (t - t)) = exp [Ft] ■ B. 



(25) 



Around the steady-state, the process is stationary, which 
means the correlation function depends upon time dif- 
ference only. Also note that the characteristic corre- 
lation time t c = ||r|| _1 is related to the Jacobian T 
of the deterministic equations, and therefore cannot be 
divorced from the deterministic relaxation time. As a 
consequence, representing the fluctuations a(t) as white 
noise (r c — > 0) is not justified. 

The great advantage of the linear noise approximation 
is that the autocorrelation function of the steady-state 
fluctuations can be calculated directly from the macro- 
scopic reaction rates in an algorithmic fashion 21]. Fur- 
thermore, since T and D arc derived from the known 
propensity and stoichiometry of the reactions, the statis- 
tics of a are fully determined and are not tunable by 
some ad hoc prescription. 



V. MEAN FIRST PASSAGE TIME 

Bistability is a property exhibited by deterministic sys- 
tems. In a stochastic context, bistability is sometimes 
assigned to an equilibrium probability distribution with 
two maxima, irrespective of their separation. A more 
practical criterion for bistability is that the two states 
are long-lived and that the mean escape time from one 
state to the other is longer than the natural timescales in 
the problem. For the single- variable autoactivator model, 
we are able to compute the escape time by an explicit 
(though approximate) expression (see [31| or p. 139 
of I30j for details). Under fairly unrestrictive assump- 
tions [541 ] , the Master equation may be approximated by 
the nonlinear Fokker-Planck equation, 



dP(a,t) 







1 d 2 



= - T(a)P(a,t) + -—D(a)P(a,t), 
at oa 2 oa z 

where the functions T and D are the nonlinear analogues 
of the coefficient matrices T and D generated by the lin- 
ear noise approximation shown in the previous section. 
For our autoactivator example, the coefficients are given 
by, 

r(a) = 7 • g(a) — S ■ a -D(a) = 7 • b ■ g(a) + S ■ a. 



The nonlinear Fokker-Planck equation has no general so- 
lution for systems of dimension greater than 1, and even 
the stationary solution is often impossible to calculate ex- 
actly for such systems [55[ . In the reduced autoactivator 
model, we are fortunate to have a system with one inde- 
pendent variable, so we can write the stationary solution 
of the Fokker-Planck equation explicitly as, 



N 
D(a) 



exp 



r(q') 

D(a'Y 



da' 



where Af is the constant of normalization (see p. 124 
of [30|). Furthermore, we can explicitly write the first 
passage time r from the HIGH state to the L OW state 
or vice- versa. 



1>(x)J D(y) 



dydx 



a m id 



^{x)J D{y) 





dydx, 



where a mi( i is the unstable equilibrium point separating 
the HIGH and L OW states a* HI and a* LO: respectively. 
The function ip(x) is given by, 



ip(x) = exp 



r(*') 

D{x') 



dx' 



(see p. 139 of [30| for additional details). 

In the main text, we discuss m.m[TLO--,Hi, t hi^lo] 
along the stability curves predicted by the effective eigen- 
values. For A b = 0.1, mm[T L0 ^Hi,T H i^Lo] =8 ± 4, 
where time has been scaled to protein lifetime (^ _1 ). 
For Ah = 0.2 and Ab = 0.3, mm[TLo^Hi ,thi^lo] = 
5.6 ± 1.4 and 5.9 ± 0.3, respectively. 



VI. DETAILS OF GENETIC CIRCUIT 
EXAMPLES 

A. The autoactivator 

We describe the transcription of the activator mRNA, 
m a and the translation of activator protein A as two dif- 
ferential equations using the activation function g to de- 
scribe the time-averaged state of the promoter, 



dm a 
dt 



7m • g{A) - 5„ 



dA 



6 P A. 



(26) 



Here 7 m is the transcription rate, 7 P is the translation 
rate, 5 m and 5 P are the rates of mRNA degradation and 
protein degradation, respectively. We make the assump- 
tion that the mRNA turnover is much faster than the 
timescale of protein degradation (i.e. S m 3> S p ). In that 
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way, we justify setting the mRNA concentration to its 
equilibrium level, 



m*(A) = Y 9(A), 

reducing the model to a single equation, 
dA _ 7 m • 7 P 



dt 



g{A) - S p A, 



(27) 



(28) 



at the expense of lumping transcription and translation 
together. Re-writing the constants 7 = 7 ^" 7p and S p — S, 
we are left with the evolution equation as written in the 
main text, 



^ = 7 • g(A) -6- A, 



(29) 



where 7 is the fully activated rate of protein synthesis 
and S is the rate of protein degradation. 



1. Transcriptional activation 

The lumping together of transcription and translation 
comes at the expense of obscuring translational ampli- 
fication of the mRNA. The translational burst size is 
approximately equal to the averaged number of protein 
molecules synthesized during the lifetime of the mRNA, 
b = j 2 - @> Hfij], so we see the production term in the 
macroscopic equation is actually (bx transcription rate), 



dA , 



9(A) -§■ A. 



(30) 



In the deterministic model, the distinction between re- 
action rate and reaction stoichiometry is immaterial, but 
that is no longer true when we calculate the intrinsic fluc- 
tuations. Writing the production and degradation stoi- 
chiometry explicitly as in the main text, 



bursty synthesis: A — U A + b; 
linear degradation: A -^-> A — 1 ; 



. -9(A), 
v 2 =S -A, 



(31) 



leading to the propensity vector u = • g(A), S ■ A] and 
stoichiometry matrix S = [6,-1]. We can easily calculate 
the coefficient matrices T and D, 



[7 • g'(A) - 5} 



T> = [b- 1 -g(A)+S-A]. (32) 



It is a simple task to then determine the steady-state 
correlations of the fluctuations, 



ID 1 [b ■ 7 • g(A*) + S ■ A*) 
2~r ~~2 [7 • g'(A*) - 6} 



(33) 



which is positive since the deterministic eigenvalue 
A = [7 g'(A*) — S] < in the stable regime where the 
analysis is carried out. We write the fractional devia- 
tion 77 of the steady-state fluctuations in A as, 



V(A 2 
A* 



(b + 1) 



1 



2[1 -A g' (A*)]V A -V ce u-g(A*y 



where A* is the steady-state activator concentration and 
Aq = % is the fully-activated protein concentration and 
lu~ 2 = Vceii is the cell volume. Provided the HIGH and 
L OW equilibrium points are well-separated (g' (A*) w 0), 
we can write, 



Vlo = 



(b + 1) 



f 



A ■ V C( 



mi 



V 7 /, (34) 



where / is the fold activation. Not surprisingly, the rela- 
tive fluctuations around the LOW state are large since in 
that state, the molecule numbers are small. More impor- 
tantly for the present discussion, we see that the mag- 
nitude of the relative fluctuations depends directly upon 
the burstiness b. To determine the effect of the bursti- 
ness upon the averaged stability, we calculate the stabil- 
ity matrices j' ^ and (where time has been scaled 
with respect to the protein lifetime: t —* t ■ S^ 1 ), 

J(°) - [Aq g' A (a) - 1] W J« = [A g'^a)]w a(t), 

from which the Laplace transform of the autocorrelation 
function J c (s) is derived, 

00 

w 2 J c (s) = lu 2 [A g"f [ (a(t)a(0))e[ Ao9 '- 1 ]*e- st di. 



Referring to Eq. [25l the steady-state fluctuations have 
exponential time-autocorrelation function so that the in- 
tegrand becomes, 



2i , \ 2 m "1 2 ( b+1 ) A o9 

uj J c (s) = -uj [A g \ 



2 [A g> - 1] 

{A„g'-l}t e {A g'-l}t e -st dt 



(35) 





Evaluating the integral, 
2 * , v (b + 1) Alg[A g'fK A 



uj J c (s) = -- 



1 



2 K A [A g> - 1] A s - 2 [A g' - 1] 

(36) 



From the stability matrices, we are able to calculate the 
approximation of the effective eigenvalue A' from Eq. I19| 



A' = [A g' - 1] + 



co 2 (b + l) K A A^jg'Tg 
K A 2 [A g' - l} 2 



(37) 



where we identify ui~ 2 = V ce u as the volume of the cell. 
Collecting the constants into groups, we write the the 
effective eigenvalue A'(^4*) as, 

A' = A + ^-Korr = A jl - A 6 • h g (A*] 

(38) 

where Ah = ^t, 1 ^ w \ r — is the discrete change in reac- 
tant molecule numbers, scaled with respect to the num- 
ber of activators required to initiate activation (K A ■ 
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Vceii), representing the relative change in protein num- 
bers incurred by the stochastic reaction events. (In a 
sense, Ka represents the characteristic concentration of 
the activator: for activator concentrations far less than 
Ka, there is no activation and for concentrations far 
above Ka, the promoter is fully activated.) The second 

term in Eq. EH h ^,g(A*)^j = ^ A °\(^) 9 contains 

the details of the regulatory mechanism [4j and depends 
strongly upon the stability of the deterministic system 
through A. It is the interplay between the fluctuations 
(through A&) and the macroscopic stability of the steady- 
state (through h) that ultimately decides the averaged 
stability of the stochastic system. 



2. Accuracy of ESA 

To compute the accuracy of the effective stability 
approximation as a function of the molecule numbers 
for the translational autoactivator model, the corrected 
eigenvalue A' computed above (Eq. [38J) is compared 
to the short-time Lyapunov exponent of the ensemble- 
averaged perturbation modes computed by stochastic 
simulation 

For a system slightly perturbed from the steady-state 
x s , the short-time Lyapunov exponent (A) is defined as, 

lim In \ (x p (t)) — x s \ = const. + (A) • t. 

A numerical calculation of (A) is obtained by taking 
the ensemble average (over an ensemble of 10 5 mem- 
bers) of x p (t) determined by stochastic simulation. The 
slope of the natural-log difference between the numeri- 
cally generated perturbation mode and the steady state, 
hi\(xp(t)) — x g \, is fit by linear regression over a time 
span corresponding to the protein lifetime (i.e. S^ 1 = 30 
minutes). To compare the stochastic simulation with 
the ESA, we focus upon three points in the parameter 
space of the autoactivator (Figure la, filled circles) - one 
point well inside the bistable regime (-^A- = 2.5, / = 
80; red), one near the boundary predicted by the ESA 
{■g^ = 3.5,/ = 80; green), and one well inside the 

monostable regime (-^ = 5, / = 80; blue). Figure lb 
compares the resulting Lyapunov exponent (A) (dashed 
lines) with the ESA prediction A' (solid lines), where the 
line colors correspond to the colors of the filled circles 
in Figure la. Here, the burstiness in protein synthesis 
is held constant at b = 9, and the characteristic num- 
ber of molecules in the system, Ka • V ce u, is increased 
from 5 to 50. (In the main text, Ka • V ce u — 25 so 
that a burstiness of b = 9 gives a discreteness param- 
eter of Ab A = Ka-v a — 0-2.) As the number of 
molecules in the system is increased, the ESA and the 
numerical simulation results converge. The figure shows 
the effective stability of the transcriptional autoactivator 
model is well-characterized by the ESA for systems with 
K A ■ Vcm > 20. 



A 

K A 




10 20 30 40 50 60 80 100 

Fold Activation / 




10 20 30 40 50 
K A' V cdi (molecules) 

FIG. 4: Accuracy of the effective stability approximation 
(ESA) as a function of the number of molecules. (A) Focusing 
upon three points in the parameter space of the autoactivator 
model (see Figure 2a in the main text), it is possible to com- 
pare the ESA with the results of numerical simulation. (B) 
The short-time Lyapunov exponent of an ensemble average of 
the perturbation modes about the LOW state (dashed lines) 
approach those values of A' predicted according to Eq. 25 
(solid lines) for systems with increasing values of Ka ■ Vcell, 
which specifies the order of molecule numbers to turn on/off 
the gene. Here, the burstiness of protein synthesis is held 
constant at b = 9, and each data point is computed from a 
sample of 10 trajectories - colors of the curves correspond to 
the filled circles in panel A. 



3. Translational activation 



To model the translational activity, we redefine the 
transcription rate to be constant j, where b is the maxi- 
mum burst size at full activation, and allow the activator 
to control the translation rate through the stoichiome- 
tery. We write the synthesis and degradation reactions - 
in analogy with Eq. 1311 above - as, 



.4 



.4 



A + b ■ g (A) 



V2 - 



— 1 
~ b ' 

- S ■ A, 



(39) 
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where the translational activation affects the stoichiom- 
etry through the synthesis step-size b ■ g(A). Notice that 
the deterministic equation &rr — S ■ i> = Aq g(A) — A is 
identical to the deterministic equation for the transcrip- 
tional autoactivator in the previous section. Nonetheless, 
the change in synthesis stoichiomctry from b <—> b ■ g (A) 
has a noticeable effect on the resulting stability. As 
above, we calculate the effective eigenvalue, 



A' = A 1 - 



(b-g(A*) + l) 



1 



V a 



K A 

where h( ■ ) is as in Eq. [38J The difference from the 
transcriptional case is that the burst-size itself is atten- 
uated in the LOW state, and the discreteness parame- 
ter approaches the minimal value A& — > l/(2V ce u ■ Ka), 
thereby increasing the residence time in the LOW state. 



Identification of dimensionless parameters in the deter- 
ministic model comes from considering the rate equa- 
tions, 



d 
dt 



A 
R 

C 



= Su = 



(42) 



1A ■ 9 



A 

K A ' 



ir ■ g 



Sr-R- 
AR- 



■ A — k c ■ A- R 

- k c ■ A- R + 5 A 
Sa-C 



C 



B. Genetic oscillator 



The parameters of Vilar et al. [12[ correspond to the 
reduced model parameters: 

7A = 25 nM h- x ,K A = 0.5 nM, f A = 10, (40) 

lR = 5 nM h-\K R = 1 nM, f~ l = 0, 

K C = 2 x 10 2 nM' 1 h~ l , and S A = 1 hT 1 , 

where, for simplicity, we make the approximation that 1 
molecule / 1/im 3 « 1 nM and set V ce ii = 100/im 3 . Fur- 
thermore, the mRNA degradation and translation rates 
in the original model give an activator burst size of bA = 5 
and a repressor burst size of b R = 10. 

1. Details of the stochastic model 

The reduced model (Eq. 6 in the main text) is com- 
posed of six elementary reactions: 



b A 



■ 9 



A^ A + b A 

A -> A - 1 v 2 =8 A - A 

{A,R,C) -> (A- l,R- 1,C + 1) u 3 = k c -A-R 

R^R + bR Vi = ft. g (£. )f 

R^R-1 v 5 =5 R -R 

{R,C)^(R+1,C-1) v & = 5 A -C 
The stoichiometry matrix S and the propensity vector u 
arc then written as, 



S = 



b A -1 -1 
-1 b R -1 1 
1 -1 

6a- A 
k c ■ A - R 

9 (A./*) 

Sr-R 
Sa-C 



(41) 



In what follows, it will be convenient to call 7 

— 7A 



and 



1A 

A = Scaling the concentrations with respect to the 



characteristic concentration Aq (i.e. A = A' ■ Aq, etc.) 
and time with respect to the activator lifetime, t = t' -Sa, 
the rate equations become, 



1-9 



d 


A' 




R' 


dp 


a 


)- 


A 1 - 


Sr 
Sa 


■R' 



K-C-Aq 

Sa 



Sa 



Sa 



A'-R'- C 



■A'-R' 
■ A'-R' 



(43) 



a 



The two additional dimensionless constants are the scaled 
rate of dimerization k — h ' c s '^° and the ratio of the re- 
pressor and activator degradation rates e = |^-. Hence- 
forth, the primes denoting the dimensionless quantities 
will be dropped. 



Since the variance in the fluctuations is found from the 
auxiliary matrices T and D (cf. Eq. 24), and T is the 
Jacobian of the deterministic system, the dimensionless 
stochastic parameters are most easily found by consider- 
ing D = S • diag[i/] • S T , 
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b A ■ 1A ■ 9 a + 5 A ■ A + 7c • A ■ C 7c ■ A ■ C - 7c A-C 

D = -ic-A-C b R - lR - g R + 8 R - R + lc ■ A-C + 8 A -C - 7c • A ■ C - S A ■ C 

--ye -A-C -7 • A-C -S A -C lc -A-C + 5 A -C 



where gt = g (-j^Tifij- As above, we scale the concen- 
trations with respect to A and divide through by 5 A . 
Evaluating D at the steady-state (A* , R* ,C*), where 
M = = — Q 5 provides the additional simplifi- 
cations derived from the rate equations above, written in 
dimcnsionlcss form, 



g A = A* + k ■ A* ■ R*, 
1-9R + C* = e • R* + k ■ A* ■ R*, 
C* = k ■ A* ■ R*. 



(44) 



Hence, the matrix D is written in terms of reactant num- 
bers as, 



D 



7 • A 



2 

c* 
-c* 



a a 

1 



C* 



-c* 



(b R + l) 



g R + 2C* -2C* 
2C* 2C* 



(45) 



Comparing each diagonal element with the characteristic 
mean reactant number of that species (N A ~ K A V ce u, 
N R ~ K R Vceii), and ignoring parameters coming from 
the deterministic model {g A , g R , and 7), we have three 
additional constants - the discreteness in the activator 
number A- - 



! , the discreteness in the 
! and the extent 



2 K A -V € 

repressor number Af,^ = ( - b " R 2 +1 j^—j- 

of dimerization K ^. v — ■ In the main text, we focus upon 
the effect of varying trie deterministic parameter e and 
the stochastic parameter A{, A . 



VII. ALGORITHMIC IMPLEMENTATION OF 
THE THE EFFECTIVE STABILITY 
APPROXIMATION 

The corrections to the deterministic eigenvalues are 
computed by solving the resolvent equation for the the 
effective eigenvalues A', 



det[A' -I- J<°) 



1 



-Jc(A')], 



(46) 



(Eq. 12 in the main text). In this section, we provide 
a step- by-step algorithm to form the matrices J^ ^ and 
J C (A') from the deterministic reaction rates. In the fol- 
lowing, the deterministic state vector is denoted by x and 
ol denotes the fluctuations in each of the components of 
x (c.f. Section I-C above). The first three steps of the al- 
gorithm come from the paper by Elf and Ehrenberg [2l| . 

1. Write the various reactions in terms of their propen- 
sity and stoichiometry. The deterministic reaction 
rates are formed by the product S • u {cf. Eqs. 31 
and 41 above). 



2. From S and u, construct the matrices T and D, 
d[S ■ u]i 



r«(x) 



D(x) = S • diag[i/] • S J 



(47) 



3. Compute the steady-state covariance in the fluc- 
tuations a. by solving the fluctuation-dissipation 
relation for each of the entries in the symmetric 
covariance matrix S (where 3y = Hj*, = {en Oj}), 



r(x s )-H + H-r J (x s ) + D(x s 



0. 



(48) 



The steady-states x s are calculated from the de- 
terministic reaction rates by solving the algebraic 
equations ([S • i/] x=Xs ) = 0. 

Evaluated at the steady-state, the fluctuation- 
dissipation relation is simply a ^d{d + 1) system of 
linear equations that determine the symmetric en- 
tries of H (where d is the dimension of the system) . 
For more details regarding the general solution of 
the fluctuation-dissipation relation, see [l8l |. 

4. Compute the matrices 3^ and jW(t), 

j(°> = r(x s) j(i) (f) = ^ + ^W) u (49) 

ouj 

5. Calculate the matrix J c (i), 

J c (t) = (j( 1 )(t)-exp[jWt]-J«(0)), (50) 

where exp[j(°)i] is the matrix exponential of 
The matrix J c (£) will be composed of lin- 
ear combinations of the autocorrelation functions 

th 



(ai(t) cej(Q)). Replace each of these by the 
element of the matrix exp[j(°) t] ■ H, 

(oi(t) a,(0)) = [exp[jW f].E] 
{cf. Eq. 25 above). 



n > 



(51) 



6. The correction matrix 3 c {t) is composed of expo- 
nential terms of the form e°', facilitating the com- 
putation of the Laplace transform J C (A'). Simply 
replace each term e at with (A' — a) -1 , 



Jc(A') = J c (*)| e at_( A /_ a )-l 

7. Solve the resolvent equation for A', 

1 



det[A' ■ I - J(°> 



V, 



cell 



-Jc(A')]. 



(52) 



(53) 



The algorithm described above is easily implemented 
in symbolic mathematics packages. A version coded in 
Mathematica is available from the authors upon request. 
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